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DIFFUSION MODEL STUDY IN CHEMICALLY REACTING AIR 
COUETTE FLOW WITH HYDROGEN INJECTION 1 

By Randolph A. Graves, Jr. 

Langley Research Center 

SUMMARY 

An analytical study of the effects of hydrogen injection and chemical reaction on the 
flow properties of Couette flow has been conducted. Special emphasis was given to the 
diffusion model assumed for the calculations. Three diffusion models were chosen for 
the analysis: Fick's law (binary diffusion), multicomponent diffusion, and an approxi- 

mation to the multicomponent diffusion. In the Fick's law model, three methods of 
obtaining the diffusion coefficient were also investigated. 

Implicit finite-difference numerical solutions to the governing equations for Couette 
flow were obtained for the three diffusion models over a range of hydrogen injection 
rates. The results indicate that there are significant differences between the solutions 
for the diffusion models and these differences are manifested most in the concentration 
profiles and the wall heating rates. 


INTRODUCTION 

The use of mass-transfer cooling to reduce aerodynamic heating encountered in 
reentry thermal environments has become widely accepted. Whether this mass-transfer 
cooling is accomplished by ablation or transpiration, the gases injected into the boundary 
layer are generally very different from those in the main stream flow. Since the con- 
vective heating reduction, as shown in reference 1, is greatest with low-molecular-weight 
gases, molecular hydrogen is usually a major component of the injected gases especially 
in the ablation of polymeric materials. The introduction of hydrogen into boundary-layer 
flow complicates the analysis because large property variations occur and molecular dif- 
fusion and chemical reactions must be considered. 


Ipart of the information presented herein was included in a thesis entitled 
"Chemically Reacting Couette Flow With Hydrogen Injection for Two Diffusion Models" 
submitted in partial fulfillment of the requirements for the degree of Master of Science 
in Mechanical Engineering, Virginia Polytechnic Institute, Blacksburg, Virginia, June 
1969. 



In most analyses, Fick's law (binary) diffusion (ref. 2) is assumed since it is a sim- 
ple and easily applied approximation to the exact (thermal diffusion being neglected) but 
mathematically cumbersome Stefan -Maxwell multicomponent diffusion model (ref. 3). 
Recently, a third diffusion model, which is a more accurate approximation to the multi- 
component diffusion model, has been proposed (ref. 3); this model utilizes a bifurcation 
of the binary diffusion coefficients to allow explicit solution of the Stefan- Maxwell rela- 
tions for the diffusive fluxes. However, since both the Fick's law and bifurcation models 
are approximations, the calculated diffusion velocities may be in error, especially when 
there are large differences in the molecular weights of the diffusing species as is the 
case when hydrogen is present in an airstream. Thus, a comparison of the diffusion 
models is necessary to provide an estimate of the errors incurred in using the approxi- 
mate models when low-molecular-weight gases diffuse through heavier gases. 

There exists little information in the literature concerning the effects of the diffu- 
sion model on the solutions obtained for a chemically reacting airflow with hydrogen 
injection. There are no direct comparisons between the approximate diffusion models 
and the exact multicomponent diffusion model available from the literature. The analysis 
of Libby and Pierucci (ref. 4) does consider hydrogen injection into a laminar air bound- 
ary layer with variable properties, a chemical reaction, and multicomponent diffusion, 
but these solutions are compared with rather limited (Prandtl and Schmidt numbers equal 
to 1) solutions and give no insight into the effect of the diffusion model utilized. The 
present analysis differs from the analysis of reference 4 in that the approximate diffusion 
models employ the same assumptions as the multicomponent diffusion analysis, except for 
the diffusion model itself. 

In making a comparison of the diffusion models, any simplification that can be used 
without concealing the important aspects of hydrogen injection into an air boundary layer 
is desirable. In the literature the one-dimensional Couette flow has been used as a 
simulation of the two-dimensional laminar boundary layer (refs. 5 and 6); however, the 
sources available consider only hydrogen injection into an air Couette flow with constant 
properties and no chemical reactions. The principal analysis is that of Eckert and 
Schneider (ref. 5), but because of their assumptions of no chemical reactions and incom- 
pressible Couette flow, their solutions are of limited usefulness. A variable property 
analysis is given in reference 6 where hydrogen is injected into a nitrogen stream, again 
with no chemical reactions and for binary diffusion only. 

The present analysis differs from those of references 5 and 6 in that variable 
properties, a chemical reaction, and three diffusion models are considered. Also, the 
present analysis does not employ the flame-sheet approximation as did Libby and 
Pierucci to define combustion but instead a diffusion flame results from the solution of 
the governing equations. 
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The primary purpose of this study is to compare the results obtained from the use 
of the three diffusion models. Also, as a result of this study, the effects of variable 
transport and thermodynamic properties and a chemical reaction on Couette flow can be 
observed. Finally, the influence of various methods of evaluating the Fick's law diffu- 
sion coefficient can be observed. 

As in the references cited, the one-dimensional Couette flow model is used as an 
approximation of the two-dimensional laminar boundary layer; however, it is recognized 
that under the conditions of the present analysis, this approximation is not accurate, but 
the Couette flow model does allow a vehicle by which the diffusion models can be com- 
pared. In this Couette flow representation, the velocity of the moving plate represents the 
free-stream velocity, whereas the distance between the plates simulates the boundary- 
layer thickness. 


SYMBOLS 


Cp specific heat of gas mixture 

Cp j specific heat of individual species 

© binary diffusion coefficient 

D Fick's law diffusion coefficient 

D average diffusion coefficient 

F w nondimensional shear stress at wall 

f diffusion factor 

h static enthalpy 

K mass fraction 

Kp equilibrium constant 

M molecular weight 

M m mixture molecular weight 
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N 

NMa 

N t 

P 

Qw 

q w 

R 

Re,, 


T 

U 

u 


X 


a 


finite-difference station number 


Mach number 


total number of finite -difference stations 


pressure 


heat -transfer rate into wall 


nondimensional heat -transfer rate into wall 


universal gas constant 

injection parameter of reference 6, pvs f — dn 

J 0 M 

distance between porous surfaces 


temperature 

nondimensional temperature for Lennard-Jones collision integral 

dimensionless flow velocity 

flow velocity 

diffusion velocity 

mass average velocity 

mole fraction 

coordinate normal to lower porous surface 
pseudo mass fraction (eq. (42d)) 

p 7 ? dr? 

coordinate parameter of reference 6, \ — - 

J 0 V- 
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reference coordinate of reference 6 


■ i 1 


dr] 


0 M 


coefficient used in bifurcation diffusion model calculations (eq. 
coefficient used in bifurcation diffusion model calculations (eq. 


nondimensional coordinate of reference 6, -5L 

’ «€ 

nondimensional mass addition rate 


maximum energy of attraction 


number of atoms of element k in a molecule of species i 


nondimensional coordinate 


nondimensional temperature 

T - T 

nondimensional temperature of reference 6, w 


Tco " T w 


total thermal conductivity 


translational thermal conductivity 


internal thermal conductivity 


viscosity 


number of species (4) 


mass density 


collision diameter 


shear stress 


reduced collision integral for diffusion 


(42b)) 

(42c)) 
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CO 


species production rate 


<p coefficient used in viscosity calculation (eq. (26)) 

i// coefficient used in thermal conductivity calculation (eq. (29)) 

Subscripts: 

i,j ith or jth species 

k kth element 

m gas mixture 

o no injection 

w wall (lower porous surface) 

free stream (upper porous surface) 
ref reference condition 

A tilde ~ over a symbol denotes an elemental. 

ANALYSIS 

Figure 1(a) shows the one -dimensional Couette flow model used in the present 
analysis and figure 1(b) gives the corresponding finite-difference representation. The 
lower porous surface, at y = 0, is stationary whereas the upper porous surface, at 
y = s, moves with a uniform velocity u^. The lower surface is at the temperature T w 
and the upper surface at Too. The hydrogen gas, initially at temperature T w , is injected 
uniformly and perpendicularly into the flow through the stationary surface, and is removed 
uniformly through the upper surface in concept only since the boundary conditions require 
that the hydrogen concentration be zero at the upper surface. 

Equations of Motion for Couette Flow 

By use of the assumptions of reference 6, the basic governing equations of motion 
for Couette flow can be reduced to the following forms: 
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Continuity: 


d(pv) _ 0 
dy 


Momentum: 


pv 


du 

dy 


dy\, dy ) 


Energy: 


pv — = — (x Y p.V jhi 

dy dy \ dy ) l^dy J dy L 1 1 


Species continuity: 

diq 

dy ' dy 


dK* 

pv *r + Jj( p i v 0 = "i 


(1) 

( 2 ) 


(3) 


(4) 


A simplification of the species continuity equations can be obtained through the 
introduction of the concept of elemental mass fractions as expressed by Lees in refer- 
ence 7. The elemental mass fraction concept results from the fact that the mass of 
individual chemical elements is preserved in any chemical reaction not involving nuclear 
transformation. The elemental mass fraction is given by the expression: 


K, 


M k 


= £ 5 *m? k ‘ 


The elemental continuity equations can be obtained by multiplying equations (4) by 
and summing over i, and, as a result, the elemental equations 


(5) 




M k 


ik Mi 


dK k d 

PV -dT + dy 



M 




(6) 


are obtained. The introduction of the elemental mass fraction eliminates the species 
production terms aq of equations (4) and reduces the number of calculations to be made. 
There is now one equation of this form for each element as opposed to one equation of the 
form of equations (4) for each chemical species. 


In the present analysis there will be three elements H, N, and O, and four chemical 
species O 2 , H 2 , N 2 , and H 2 O considered with one chemical reaction of the form: 


h 2 + |°2 


h 2 o 
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This same chemical system was used by Libby and Pierucci in reference 4, and does 
not consider dissociation or ionization. The maximum gas temperature in the present 
study is less than 2400° K, and at this temperature and the pressure of 1 atmosphere 
assumed for this study, the amount of dissociation of O 2 , H 2 , and N 2 is negligible 
(1 atmosphere = 101.325 kN/m2). The species considered have the necessary variation 
in molecular weight which is essential to the diffusion-model comparisons. 


Boundary Conditions 

At the moving surface (y = s), the following boundary conditions apply: 


T = T 


OO 


At the wall (y = 0), the boundary conditions are: 
u = 0 
T = T W 

The boundary conditions on the elemental mass fractions are derived as follows. Inte- 
gration of the continuity equation (eq. (1)) yields 

pv = Constant = (pv) w 

By using this relation, the elemental continuity equations can be integrated to give 


pvK k + 2/ ik “ p i V i = Constant (?) 

The following subscript notation is adopted for the elements: 


Element 

Subscript 

O 

k = 1 

H 

k = 2 

N 

k = 3 


By considering the injected hydrogen first, equations (7) become: 
P 2 V + PgV 2 = Constant 

Pg(v + V 2 ) = Constant 

p 2^2 = Constant 
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Evaluating the constant at the wall (y = 0) yields: 


= (?2^) w = (P V K 

Thus, the wall boundary condition on the elemental hydrogen mass fraction becomes 


or 


(pv)wK k ,w + (y < ik ^ - <pv) 


(k = 2) 


(k= 2) 


(8) 


A similar procedure is followed in evaluating the constant for the main stream 
elements where 


5 1*1 m {h V l) v = ° 

V3 = - 0 

The boundary conditions for these elements are: 


(k= 1, 3) 


(9) 


In order to simulate the two-dimensional boundary layer, the elemental mass fraction for 
hydrogen must approach zero at the upper boundary. This condition creates a corre- 
spondingly small elemental hydrogen density and since the elemental continuity equation 
pv = Constant must be satisfied, the elemental transverse velocity becomes very large. 
This condition also introduces some uncertainty since the transverse velocity was 
assumed to be small in comparison with the main flow velocity to make possible the 
reduction of the general equations of motion. However, the inaccuracies incurred are 
confined to the region immediately adjacent to the upper boundary and are inherent in the 
use of the one -dimensional Couette flow to simulate the two-dimensional laminar bound- 
ary layer. 


Nondimensional Form of the Governing Equations 
The following new variables are introduced: 


V 


= 1 




i 




6 = P vs = (PV) w S 

Moo M oo 

The governing equations in nondimensional form are 
Momentum: 

6 dU _ d f li dU' 


t (10) 

dr) dr]\^oo dr, j 

Energy: (in addition to the nondimensional variables, the free-stream specific heat 

C is needed to provide the following nondimensional energy equation.) 

P>°° 


dh _ d / X d6\ u ~ M /dU\ 2 


C p>00 T M d ^ ^(Cp^oo dr, 


Cn M 

P j°o 00 


00 


dv) 



Pi ViS 


O Cp )co T w 


(ID 


Elemental continuity: 

d *k _ dL M k p i V i S \ 

dr) drjWL/ 1^ Mj Moo I 


( 12 ) 


Nondimensional Boundary Conditions 
The nondimensional boundary conditions at r, = 1 are 
U = 1 

0 = 1 

= ^k,oo 

At 77 = 0, 

U = 0 
0= 0 W 

(k = 2) (13) 

(k = 1 , 3) (14) 
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Heat Transfer Into Wall 
The heat -transfer rate into the wall is 


«w - §• - 2 PiVA 


w 


Transformation of equation (15) yields 

Qw s / x m d0 V P i ViS h i 




w McoCp^Tco ^ooCp ?00 dr/ Zv Moo Cp ?0O Tooy 


-z- 


Shear Stress at Wall 
The shear stress at the wall is 
du 

T W = 

w m dy w 

Transformation of equation (17) yields 
T w s ^m dul 


F = 

w Moo d77 


w 


(15) 


(16) 


(17) 


(18) 


Gas Properties 

The chemical thermodynamic and transport properties are calculated by the methods 
given in this section. The gas mixture is assumed to be at a pressure of one atmosphere 
for all calculations; however, comparison cases determined at lower pressures indicate 
only a minor influence of pressure on the solutions. 


Chemi cal composition .- The following reaction is considered for the present 
analysis: 

H 2 +!°2=H 2 0 


In addition, N 2 is present in the main Couette flow; thus, there are four chemical species 
to be considered in the equilibrium calculations. The equilibrium constant is related to 
the mole fractions by 


Kp-P -^ 2 


X 


h 9 o 


Xh 2 (Xo 2 ) 


172 


Substitution of 


K . _ X i M i 
1 M m 


(19) 


(20) 
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into equation (5) yields 



Mk 

M m 

i=l 


1/ 

1 Vi 


The relation 


(k=l, 2, 3) 


(21) 


l 

i=1 


x i = i 


( 22 ) 


combined with equations (19) and (21) constitutes a system of five equations in the five 
unknowns, Xi (where i = 1, 2, 3, 4) and M m . These equations were combined to 
eliminate four of the unknowns; one equation remained to be solved numerically for the 
mole fraction Xq 2 * The equilibrium constant used in these calculations is taken from 
the JANAF tables (ref. 8). The species mass fractions are determined from the mole 
fractions by equation (20). 


Thermodynamic properties .- The mixture density is obtained from the equation 
of state 


_ P M m 
p m rt 


(23) 


and the enthalpy of the individual species is taken from the JANAF tables (ref. 8) and the 
mixture enthalpy is calculated by 


‘m 


-l 


Kihi 


i=l 


(24) 


Transport properties .- Rigorous kinetic theory expressions for the viscosity and 
thermal conductivity of gas mixtures have been developed and are presented by 
Hirschfelder, Curtiss, and Bird in reference 9, but these expressions are mathematically 
cumbersome. Somewhat simpler relations, which are approximations derived from the 
rigorous expressions, are given by Brokaw in reference 10 and are used in the present 
analysis. A comparison and discussion of approximate and rigorous expressions for an 
equilibrium reacting gas can be found in reference 11. In the present analysis the pure 
species viscosity and thermal conductivities are obtained from reference 12 where they 
were calculated by using the molecular constants given in table I and the Lennard-Jones 
(6-12) potential. (See ref. 9.) 

The mixture viscosity is calculated from the pure component viscosities with the 
relation 
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( 25 ) 


M = y tti 

M m / v 

'■"•K 

The coefficients (p y were derived in reference 10 by use of rigid-sphere theory and 
are a function of the pure component viscosities and molecular weight ratios 


; \ 1 /2/mA1/4 


r / mA 1 / 2 

2 *N) 

The pure component viscosities are plotted in figure 2 where it is seen that hydrogen has 
a considerably lower viscosity than the remaining species. 

The mixture thermal conductivity is obtained from the relation: 

^m = ^m + ^m (27) 

The translational mixture conductivity is obtained from the pure component translational 
conductivities with the relation 




1=1 1 + > A 


The coefficient is obtained from the viscosity coefficient cp^ by the following 

relationship obtained from reference 10: 


C (Mi - Mi) (Mi - 0.142MA ") 

*ij = 4>ij 1 1 + 2.41 ^ M-J U V (29 

3 y L ( M i +M i) 2 Jj 

The internal mixture conductivity is obtained from the pure -component internal conduc- 
tivities with the relation: 


v . M 

£•1 £ - 

i=l 


■111111111 II II llllll II IIIIH II napnpui ■ Ml a imm in hi i ■iiiiiiiiHini 


ii ■mu 


The total thermal conductivity for each species is shown in figure 3 and as with the pure- 
component viscosities, the hydrogen is again very different from the remaining species, 
its thermal conductivity being much greater. 


Diffusion Transport 

The purpose of the present analysis is to compare solutions to the governing equa- 
tions for Couette flow by using three different diffusion models: the approximate Fick's 
law diffusion model, the exact multicomponent diffusion model, and the bifurcation model. 
As will be seen below, the exact multicomponent diffusion model entails many mathe- 
matical operations and from a numerical analysis standpoint is not as desirable as the 
simpler but approximate Fick's law model. In the Fick’s law diffusion model, three 
methods of calculating the diffusion coefficient are explored, the results indicating a 
wide variation in the solutions obtained. 

Multicomponent diffusion.- The multicomponent diffusion fluxes were calculated by 
use of the Stefan-Maxwell relation from reference 9: 



Equation (31) can be rearranged to a more convenient form: 


(31) 


(32) 



XiXj 

®U 


(33) 


Multiplying equation (33) by p/Poo an d introducing the nondimensional coordinates 
yields: 



f X i X i 

PjVjS 

Poo d77 

,4 K sei: 

j M OO 


j*i 


PjVjS y XjXj 
M oo /-/ Kj 

3=1 J 


(34) 
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Similarly, multiplying equation (32) by s/jjl^ yields 


V Pi v i s 

L Moo 

i=l 


= 0 


( 35 ) 


For the ^-species gas mixture, the diffusion fluxes p^V^sjii ^ are obtained from the 
simultaneous solution of v - 1 relations of the form of equation (34) and the relation 
given by equation (35). 


The binary diffusion coefficients used in equation (34) are calculated by use of the 
following relation from reference 9: 


_ i SMilVh 

= 0.002628 A= —L 


+ Mj 


1/2 


P ( CT ij ) 2s2 i?’ 1) 


i] 

The collision cross section cr— is obtained from the relation 


ff ir 


a i + 0 j 


(36) 


where the collision cross sections for each species are obtained from reference 12 and 
are given in table I. 

(1 1 )* 

The reduced collision integral ’ is based on the Lennard-Jones (6-12) 
potential and is taken from reference 9 where it is tabulated as a function of the nondimen- 
sional temperature T*j which is defined as 


T* = T 

« - ^/k 

The maximum energy of attraction 




in °K is obtained from 


where the maximum energy of attraction for each species is taken from reference 12 and 
is given in table I. The binary diffusion coefficients obtained from equation (36) to be 
used in equation (34) are shown in figure 4, where it is apparent that the interactions 
involving the low-molecular-weight hydrogen produce larger binary diffusion coefficients. 


Fick’s law diffusion .- The Fick's law diffusion fluxes are calculated according to 
the following relation: 


dK; 

p, v i = -pD -j-! 


(37) 
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Knuth in reference 2 states that a sufficient condition for the applicability of equation (37) 
is that the binary diffusion coefficients are equal to each other and to the Field s law dif- 
fusion coefficient. This assumption makes the Fick’s law diffusion coefficient a pseudo 
binary diffusion coefficient and in the literature Fick's law diffusion is generally referred 
to as binary diffusion because of the appearance of equation (37). The term binary dif- 
fusion is adopted here for discussion purposes. 

By multiplying equation (37) by 1/p^ and introducing the nondimensional coordi- 
nates, equation (37) becomes 

(38) 

Moo dr > 


The calculation of the Fick's law diffusion coefficient can be accomplished in a number of 
ways; however, the following three methods have been selected for the present study. 


Method 1: In the first method the diffusion coefficient is assumed to be independent 
of the molecular concentrations and is given by the self-diffusion relation from 
reference 9: 


D = 0.002628 


(t 3 /m) 1,/2 


(39) 


where a, M, and e/k (needed to calculate n(-M) ) are mixture averages as given in 
table I. Thus, D is dependent only on temperature and pressure. The diffusion coef- 
ficient calculated by equation (39) is given in figure 5. By comparison with figure 4, it is 
apparent that use of the average molecular constants causes the Fick's law diffusion coef- 
ficient to lie in the region of the heavy-molecule binary diffusion coefficients; thus, little 
of the effect of the low-molecular-weight hydrogen is provided. 


Method 2: In the second method the diffusion coefficient is allowed some dependence 
on the molecular concentrations by allowing the molecular weight in equation (39) to vary 
as the mixture molecular weight. The diffusion coefficient is given by 


D = 0.002628 


(T 3 /M m ) 1/2 

pa2fi(M)* 


(40) 


where a and e/k are given in table I. There is some inconsistency in using this pro- 
cedure since the molecular weight is allowed to vary but not the other two molecular 
constants. However, the diffusion coefficient calculated by equation (40) does provide for 
a better representation of the average diffusion coefficient as seen in figure 5 where the 
diffusion coefficient covers a wide range of values more representative of the binary dif- 
fusion coefficients seen in figure 4. The upper and lower limits on the values seen in 
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figure 5 for equation (40) were determined by assuming that the stream consisted entirely 
of hydrogen (upper limit) and of air (lower limit). 

O 

Method 3: In the third method^ the diffusion coefficient is a strong function of the 
molecular concentrations and is given by 


D = 0.002628 


(* 3 ) 1/2 

P 


2 


Xj 

(Mi) 1 / 2 ^* 1 ’ 1 )* 


(41) 


where ct^, Mi, and e^/k are the molecular constants of the ith species and Xj is the 
mole fraction. Equation (41) provides the means of allowing the Fick’s law diffusion 
coefficient for mixtures a wider variation of values than did either equation (39) or (40) , 
as seen in figure 5. The upper and lower limits for equation (41) were determined in the 
same manner as the method 2 limits. 


Bifurc ation model .- The diffusion velocities are calculated by using the following 
simplified form of the Stefan- Maxwell relation (ref. 3): 


v . _ pof @2 dZj Zq - d/3 2 \ 

1 1 h \ M m dy M m dy ) 


(42a) 


where 


H - 1 

i 

^2 = M mZr 

V I i 


Z) 


MjXj 

^2 


(42b) 


(42c) 


(4 2d) 


_ 0.002628 \P7m~T 

D = I — L — ESI (43) 

p<T ref^ 1,:L) * 

and the ^-diffusion factors f^ are obtained from a least-squares fit to the exact binary 
diffusion coefficients, to be described subsequently. Equation (42a) was obtained from 
the Stefan-Maxwell relation and equation (31) by a bifurcation of the binary diffusion 
coefficients, 

2 This method was suggested by Dennis O. Allison of the Langley Research Center. 
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(44) 


®ij= D- 
f i f 3 


where D is given by equation (43). The values for f are determined by finding the set 
of diffusion factors that gives the minimum total system residual error. Differentiating 
the total-system error relation 


R =i 1 (®«Vi • d ) 2 

i=l 3=1 

with respect to f and setting the resultant expression equal to zero yields 


(45) 


fi = 



(46) 


from which the diffusion factors can be found by iteration. The diffusion factors are 
thus obtained from a least-squares fit to the kinetic-theory binary diffusion coefficients 
(eq. (36)). It is shown in reference 13 that the diffusion factors have only a weak depen- 
dence on temperature. In the present analysis, the fj are assumed to be constant at 
the values given in table II and D has been evaluated by using the values for O 2 from 
table I for the reference constants. Table II also gives a comparison of the binary dif- 
fusion coefficients calculated by equations (36) and (44) where it is seen that equation (44) 
is a good approximation to the exact equation (eq. (36)); thus, the approximate bifurca- 
tion method should represent the exact multicomponent diffusion model fairly accurately. 

The final form of the bifurcation diffusion flux relation is found by multiplying equa- 
tion (42a) by 1/jUoo and introducing the nondimensional coordinates to get 

PiV pp (h dz i , Z 1 - K i m 

4«> d!) M m djj J 1 

This equation is from a numerical standpoint easier to evaluate than the system of equa- 
tions required by the multicomponent model (eqs. (34) and (35)). 
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Computation 


The governing equations for Couette flow with hydrogen injection can be put in more 
convenient forms for numerical solution. Equation (10) is integrated and the constant of 
integration is evaluated at rj = 0. The resulting momentum equation is 


M m dU fl 1 m du\ = gu 
M oo drj yi oo d?7 


(48) 


Similarly, equation (11) is integrated and the constant of integration is evaluated at t) = Q. 
The resulting energy equation is 


^m d6 
C p,oo^°o ^7 



C P,oo T °° 


( h m - h m,w) + £ 


P i V i s 


Me 


hj 

C P,oo T °° 


V p i V i S h i \ C 77 vt *WdU\ 2 d (49) 
4 ^=0 Cp^T^ J 0 Cp )0oToo Moo (dr? / 


The solutions to the momentum and energy equations for all diffusion models are obtained 
by an implicit finite-difference numerical technique. Briefly, this technique involves 
expressing the derivatives on the left-hand side of each equation as four-point numerical 
differences and evaluating the right-hand side at each finite -difference station. The 
resulting system of linear algebraic equations is expressed in matrix form and a solution 
obtained therefrom. 

The solution to the elemental continuity equation follows a somewhat similar pro- 
cedure. Equation (13) is integrated and the constant of integration is evaluated at r/ = 0. 
The elemental continuity equation becomes 


s _ i f ji AW 
k 5 L Mj Moo 

i=l 1 


(k = 1, 3) (50) 


Since there are three elements in the system, only two elemental continuity equations need 
to be solved since the sum of the elemental mass fractions equals unity. For the multi- 
component and bifurcation diffusion model solutions, equation (50) is solved by the method 
of successive approximations since the diffusion fluxes are given by equa- 

tions (34), (35), and (47) from which the elemental profile dependence cannot be separated. 
In the case of the Fick's law solutions, the right-hand side of equation (50) can be partially 
replaced by equation (38) and by noting the definition of elemental mass fractions (eq. (5)), 
the following elemental continuity equation results: 

6 V£i^r 0 (k=1 - 3 > <51) 
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In this simplified form, the elemental continuity equation for the binary diffusion model 
can be solved by use of the same implicit finite-difference scheme that was used to obtain 
solutions to the momentum and energy equations. The iterative solution of the finite- 
difference forms of the governing equations is accomplished by the iteration scheme 
given in reference 14. 


RESULTS AND DISCUSSION 

The results from a test case solved by the present numerical solution technique 
have been compared with the results from references 6 and 14. The flow problem is 
that of hydrogen injection into a nitrogen Couette flow for the following conditions: 

Too = 218° K; Njyj a = 12; T w = 872° K; and Re v = 0.5. In each case the solutions were 
obtained with the assumption of variable properties and exact binary diffusion. As shown 
in figures 6 and 7, good agreement was obtained between the methods. This comparison 
case contains all the essential features of the present solution technique except for the 
chemical reaction itself; hence, the present numerical technique is thought to be suffi- 
ciently accurate to carry out the present investigation. 

The remaining cases in this report consider the problem of hydrogen injection into 
air Couette flow. The values Too = 218° K; Nj^ a = 6; and T w = 872° K were held 
constant for these cases. The influence of hydrogen injection was studied by allowing 6 
to assume the values 6 = 0, 0.05, 0.1, 0.13, 0.2, 0.35, 0.5, 0.75, 1.0, and 1.3. In the 
numerical calculations the solution for 6 = 0 (no injection) was obtained with 40 finite- 
difference stations; solutions for 6 > 0 were obtained with 50 finite-difference stations. 

The no-injection temperature and velocity profiles are given in figure 8. It should 
be noted that the stream temperature increases only slightly above the wall value; thus, 
for the present conditions, the wall temperature is less than but close to the adiabatic 
wall temperature. The velocity profile is not a linear profile because of the viscosity 
variation through the stream. 

Effect of Concentration Profiles on Transport Properties 

The differences between the diffusion models are best seen in the concentration 
profiles in figure 9. It can be seen that changing 6 produces changes in both the rela- 
tive amounts of the various species and also produces variations in the profile shape. 

The biggest concentration differences occur for hydrogen, the wall concentration best 
reflecting this difference. This effect is summarized in figure 10 which gives the hydro- 
gen concentration at the wall for all the diffusion models. It is readily seen that the 
binary diffusion model concentrations are much larger than the corresponding multicom- 
ponent and bifurcation models concentrations for all values of 6 > 0.16. 
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The concentration profiles determined by the three methods of evaluating the Field s 
law diffusion coefficient show large variations in both magnitude and shape. In particular, 
the location of the flame zone is strongly dependent on the diffusion coefficient, the zone 
location moving away from the wall with increasing diffusion coefficient. 

In general, figure 9 shows that there is little detectable difference between the 
multicomponent and the bifurcation models. This result would be expected from the com- 
parison of binary diffusion coefficients given in table II. These figures also show that 
the method 3 binary model solutions, which are strongly concentration dependent, appear 
to give the best Fick's law solutions. This statement is especially true in the wall region 
which is most important for heat -transfer and shear-stress calculations. 

The concentration differences seen in figure 9 alter the mixture transport properties 
at the wall, because hydrogen has a larger thermal conductivity and lower viscosity than 
the other species. This alteration of the transport properties is readily seen in figure 11 
which gives a comparison of the mixture viscosity and thermal conductivity at the wall for 
the diffusion models. The greater hydrogen concentration of the binary diffusion model 
with methods 1 and 2 diffusion coefficients results in a mixture viscosity which is lower 
and a thermal conductivity which is higher than the other models. As would be expected 
from figure 10, there is no detectable transport-property difference between the multi- 
component and bifurcation diffusion model solutions. 

Temperatur e profiles .- The nondimensional temperature profiles are given in fig- 
ure 12. The bifurcation and multicomponent diffusion models yield essentially identical 
results for the temperature profiles over all injection rates. The results for the binary 
model solutions show only fair to poor agreement with the multicomponent profiles, and 
the method 1 binary model solutions generally give the poorest agreement. 

As an additional point of interest, the rather strong effect of the chemical reaction 
is seen by comparing the no-injection temperature profile of figure 8 with those of fig- 
ure 12. The increase in peak stream temperature over the no-injection case approaches 
a factor of three at the higher injection rates. 

Wall he ating rates .- The wall-heating-rate curves shown in figure 13 point out some 
of the largest differences between the diffusion model results. The heating rates for the 
binary model are larger than the corresponding multicomponent and bifurcation models, 
the multicomponent and bifurcation models giving essentially identical results. 

The heating rates due to conduction for the binary diffusion model with methods 1 
and 2 diffusion coefficients are generally much larger than the corresponding multicom- 
ponent solutions whereas the heating rates due to diffusion (generally negative; energy 
diffusion away from lower surface) for the multicomponent and bifurcation diffusion 
models are larger because of the greater hydrogen diffusion velocity, the net effect being 
the lower heating rates for the multicomponent and bifurcation models. 
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Velocity profiles . The nondimensional velocity profiles are given in figure 14. 
There does not appear to be any major difference between the solutions for the diffusion 
models, especially at the lower injection rates where the amount of hydrogen and water 
are substantially reduced in comparison with the oxygen and nitrogen. Again, as in the 
case of the , temperature profiles, the bifurcation and multicomponent diffusion models 
yield essentially identical results; however, the method 3 binary model solutions are also 
very close to those for the multicomponent and bifurcation models. 

Shear stress .- The shear stress for the multicomponent and bifurcation diffusion 
models is higher than the corresponding binary diffusion model solutions for all injection 
rates. (See fig. 15.) This shear-stress difference results primarily from the mixture 
viscosity variations between the diffusion models. The pure component viscosities for 
hydrogen and water are lower than those for nitrogen and oxygen; as a result, mixture 
viscosity decreases with increasing hydrogen and water concentrations as seen in fig- 
ure 11. This decreased mixture viscosity for the binary diffusion model causes the some- 
what reduced shear stress as seen in figure 15. The method 3 shear stress is not as 
close to the multicomponent solution as might be expected. 

Diffusion Coefficient Methods 

The temperature and concentration dependences of the diffusion coefficients of the 
three binary methods were investigated. A single set of concentrations and temperatures 
was provided by the multicomponent solutions at an injection rate 6 = 1.3 for which 
there is a great deal of hydrogen present in the wall region. 

In figure 16(a) it is apparent that the concentration dependence of the method 3 dif- 
fusion coefficient far outweighs its temperature dependence and also causes it to be much 
larger than the corresponding diffusion coefficients for methods 1 and 2. This larger 
diffusion coefficient in the wall region is responsible for the better comparisons of 
method 3 with the multicomponent solutions. A second multicomponent case was selected 
(6 = 0.13) in which the hydrogen concentration was much less than that for the previous 
case. The diffusion coefficients for the three methods are plotted in figure 16(b). 

Here it is seen that the method 3 diffusion coefficient is smaller than those for the 
other two methods and has the same temperature -dependent shape as these other methods. 
Figure 16(b) compared with figure 16(a) illustrates the strong effect of the hydrogen on the 
method 3 diffusion coefficients. The wider range of values given by the strong concentra- 
tion dependence of the method 3 binary diffusion coefficients is responsible for the gener- 
ally better comparisons with the multicomponent diffusion model solutions. 


22 


CONCLUDING REMARKS 


I 


A numerical study of the influence of various diffusion models on air Couette flow 
with hydrogen injection has been conducted. For this study the flow parameters were 
fixed with the moving wall temperature equal to 218° K, the moving wall Mach number 
equal to 6, and the stationary wall temperature equal to 872° K. The dimensionless injec- 
tion rate was varied between zero and 1.3. Three diffusion models were included in the 
study: the Stefan- Maxwell multicomponent diffusion model, the bifurcation model (an 
approximation to multicomponent diffusion), and the Fick's law (binary diffusion) model. 

In addition, three variations on the method of calculating the Fick's law diffusion coeffi- 
cient were investigated. 

The results of the present investigation show that the bifurcation model yielded 
essentially the same results as the multicomponent diffusion model whereas the Fick's 
law model produced results for which the agreement with the multicomponent model 
ranged from very poor to fair. The best Fick's law results were obtained with a diffu- 
sion coefficient that was strongly concentration dependent. It was also found that for low 
injection rates, the calculated heating rates and shear stresses are not influenced by the 
diffusion model, but at high injection rates, these parameters are materially influenced 
by the increased hydrogen in the airstream. 

The choice of a diffusion model for a particular problem obviously depends on the 
requirements for accuracy and ease of computation. The present case, involving hydro- 
gen injection into air, provides a severe test for the two approximate diffusion models 
studied. The results have demonstrated that when properly applied, the Fick's law and 
bifurcation models can provide a fair degree of accuracy at significant savings in numeri- 
cal complexity. It must be concluded that unless stringent requirements are placed on 
the accuracy of the results, one of the approximate models should be used. It must also 
be concluded that even when a high degree of accuracy is required, the bifurcation model 
should be investigated before resorting to the numerical complexities of the multicompon- 
ent diffusion model. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., August 11, 1970. 
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TABLE I.- MOLECULAR CONSTANTS 


Species 

e/k, 

°K 

b°<! 

. 

g/g-mole 

°2 

106.7 

3.467 

32.00 

h 2 

59.7 

2.827 

2.016 

h 2 o 

809.1 

2.641 

18.02 

n 2 

71.4 

3.798 

28.02 

Mixture 

138.503 

3.183 

20.014 


TABLE II.- COMPARISON OF BINARY DIFFUSION COEFFICIENTS 
AS COMPUTED BY THE BIFURCATION TECHNIQUE 
AND FROM KINETIC THEORY 

[t = 1000° k] 


Species 



©ij 

(eq. (36)), 
cmVsec 

fi 

©ij 

(eq. (44)), 
cmV sec 

Percent 

i 

j 

error 

°2 

H 2 

6.0174 

0.97644 

5.9018 

-1.92 

o 2 

h 2 o 

2.0200 


2.0140 

-.29 

02 

n 2 

1.5863 


1.6249 

2.44 

h 2 

h 2 o 

6.7932 

.28259 

6.9589 

2.44 

h 2 

N 2 

5.6324 

i 

5.6147 

-.31 

h 2 o 

n 2 

n 2 

1.9539 

.82811 

1.02636 

1.9160 

-1.94 




(a) Schematic diagram. 


N = Nt 


Moving boundary 


Finite -difference 



Stationary boundary 


(b) Finite-difference representation. 
Figure 1.- Couette flow model. 






Figure 4.- Binary diffusion coefficients. 
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Figure 5*- Fick*s lav diffusion coefficients. 
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Figure 7 *- Velocity and concentration profiles for 
nitrogen Couette flow with hydrogen injection. 











Binary 

model 


Method 1 

Method 2 

Method 3 

Multicomponent model 

Bifurcation model 



i 


(c) 5 = 0.35. 


Figure 9*- Continued. 
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Figure 11 Transport properties at the wall. 
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(b) 6 = 0.75. 


Figure 12.- Continued 
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(c) 6 = 0.35. 

Figure l 4 .- Continued. 
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Figure l6.- Fick ! s law diffusion coefficient profiles based on the temperature 
and concentration profiles of the multicomponent solution. 
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(b) 8 = 0 . 13 . 

Figure 16 .- Concluded 
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